Lasers passively mode-locked by means of saturable absorbers have a tendency to Q-switched mode-locking (QML). Usually, this is not the desired mode of operation and has to be avoided by careful design. This is not always easy, especially not for the low pulse energies typically associatedd with very high pulse-repetition rates.
In the following, I summarize the dynamics of such a laser and show the criterion for stable mode locking (i.e. not QML).
The following differential equation describes the dynamics of a passively mode-locked (see, e.g., C. Hönninger et al., J. Opt. Soc. Am. B/Vol. 16, No. 1, 1999):
$$\begin{align*} \dot{P} &= \frac{g - l - q_P(S)}{T_R} P \\ \dot{g} &= -\frac{g}{\tau_L} + \frac{\eta \cdot P_P}{E_{sat,L}} - \frac{P \cdot g}{E_{sat,L}} \\ q_P(S) &= \frac{\Delta R}{S} \left(1 - e^{-S} \right) \end{align*}$$,
where $S = \frac{E_P}{E_{sat,A}} = \frac{P \cdot T_R}{E_{sat,A}}$.
Let's set this up in sympy but let's - for the time being - ignore that $q_P$ is a function of $P(t)$:
In [1]:
import sympy as sym
from IPython.display import display
sym.init_printing(use_latex='mathjax')
# Variables
P, g, PP, qP = sym.symbols('P g P_P q_P', real=True)
l = sym.symbols('l', real=True)
DR = sym.Symbol('\Delta R', real=True)
EsatL, tauL, TR, eta, EsatA = sym.symbols('E_{satL} tau_L T_R eta E_{satA}',
real=True)
# Differential Equations
laserEq = [sym.Eq((g - l - qP) / TR * P),
sym.Eq(-g / tauL + eta * PP / EsatL - P * g / EsatL)]
display(*laserEq)
In [2]:
PP0 = sym.Symbol('P_{P0}')
gst_, Pst_ = sym.solve(laserEq, (g, P))[1]
Pst_ = Pst_.subs(PP, PP0)
display(sym.Eq(sym.Symbol('g_{st}'), gst_),
sym.Eq(sym.Symbol('P_{st}'), Pst_.subs(g, sym.Symbol('g_{st}'))))
To determine $q_P$, we multiply the $P_{st}$ equation by $g_{st}$ on both sides and substitute $g_{st}$:
In [3]:
qPEq = sym.Eq(P * gst_, Pst_ * g)
qPEq.subs(P, sym.Symbol('P_{st}'))
Out[3]:
Now we re-write both sides as functions of $S$ to obtain an equation for the steady-state $S$ parameter:
In [4]:
S = sym.symbols('S', real=True)
qPEq = qPEq.subs(qP, DR / S * (1 - sym.exp(-S))).subs(P, S * EsatA / TR)
sym.simplify(qPEq)
Out[4]:
Unfortunately, it turns out that this equation is rather hard to solve symbolically:
In [5]:
sym.solve(qPEq, S)
In [ ]:
dP = sym.Symbol('\delta P', real=True)
dg = sym.Symbol('\delta g', real=True)
dPP = sym.Symbol('\delta P_P', real=True)
dqPdEP = sym.Symbol("q'_P", real=True)
Pst, gst = sym.symbols('P_{st} g_{st}', real=True)
deviationEq = ['', '']
deviationEq[0] = sym.Eq(sym.Symbol('\delta Pdot'), laserEq[0].lhs)
deviationEq[1] = sym.Eq(sym.Symbol('\delta gdot'), laserEq[1].lhs)
deviationEq = [eq.subs([(P, Pst + dP),
(g, gst + dg),
(PP, PP0 + dPP),
(qP, qP + dP * TR * dqPdEP)]) for eq in deviationEq]
display(*sym.simplify(deviationEq))
Next, we drop all higher order term. Unfortunately, sympy does not (yet?) support multivariate series expansion. Therefore, we solve the problem by introducing a helper variable $\delta x$ which is subsequently dropped.
In [ ]:
dx = sym.Symbol('\delta x', real=True)
lin = [eq.rhs for eq in deviationEq]
lin = [expr.subs([(dg, dg * dx), (dP, dP * dx)]) for expr in lin]
lin = [expr.series(dx, x0=0, n=2).removeO().subs(dx, 1) for expr in lin]
lin = [sym.collect(sym.collect(eq, dP), dg) for eq in lin]
linEq = [sym.Eq(sym.Symbol('\delta Pdot'), lin[0]),
sym.Eq(sym.Symbol('\delta gdot'), lin[1])]
linEq
There we are. However, this does not look linear. There seem to be unexpected constant terms. It turns out that these are zero due to the steady-state condition:
In [ ]:
# Replace steady-state symbols with their corresponding expressions:
linEq[0] = sym.simplify(linEq[0].subs(gst, gst_))
linEq[1] = linEq[1].subs([(Pst, Pst_), (gst, gst_)])
# Rewrite in steady-state symbols for clarity:
linEq[1] = linEq[1].subs(PP0, sym.solve(sym.Eq(Pst, Pst_), PP0)[0])
linEq[1] = linEq[1].subs(l, sym.solve(sym.Eq(gst, gst_), l)[0])
linEq[1] = sym.simplify(linEq[1])
linEq
Now, we can rewrite this in handy matrix-notation in state-space form as follows:
In [ ]:
A = sym.zeros(2, 2)
for i in range(2):
for j in range(2):
col = [dP, dg][j]
A[i, j] = sym.collect(linEq[i].rhs.expand(), col, evaluate=False)[col]
B = sym.zeros(2, 1)
for i in range(2):
try:
B[i, 0] = sym.collect(linEq[i].rhs.expand(), dPP, evaluate=False)[dPP]
except KeyError:
B[i, 0] = 0
x = sym.MatrixSymbol('x', 2, 1)
u = sym.MatrixSymbol('u', 1, 1)
stateEq = sym.Eq(sym.MatrixSymbol('xdot', 2, 1), A * x + B * u)
stateEq
where $\vec{x} = \left[ \delta P, \delta g \right]^T$ and $\vec u = [\delta P_P]$. And to be complete, we add the output equation:
In [ ]:
y = sym.MatrixSymbol('y', 1, 1)
Toc = sym.Symbol('T_{oc}', real=True)
C = sym.Matrix([[Toc, 0]])
outputEq = sym.Eq(y, C*x)
outputEq
where $\vec y$ is $\left[ P_{out} \right]$.
It is now straightforward to discuss the stability of the system at its steady-state. The system is stable if the eigenvalue of $A$ all are in the left half of the complex plane. The eigenvalues of a 2 x 2 matrix
$$ \left( \begin{array}{ccc} \alpha & \beta \\ -\gamma & -\epsilon \end{array} \right) $$are:
In [ ]:
alpha, beta, gamma, epsilon = sym.symbols('alpha beta gamma epsilon', real=True)
M = sym.Matrix([[alpha, beta], [-gamma, -epsilon]])
M.eigenvals()
In the case of solid-state lasers, the radicand is negative, leading to relaxation oscillations. If $\alpha - \epsilon < 0$, these oscillations are damped, and the laser is stably mode locking. This conditions translates to
In [ ]:
A[0, 0] + A[1, 1] < 0
That's it for now. As soon as I understand limits and assumptions in sympy I may add the useful simplification for strongly saturated absorbers. The result will be:
$$ E_P^2 > E_{satL} E_{satA} \Delta R $$